A model of fasciculation and sorting in mixed populations of axons 
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We extend a recently proposed model (Chaudhuri et ai, EPL 87, 20003 (2009)) aiming to describe 
the formation of fascicles of axons during neural development. The growing axons are represented 
as paths of interacting directed random walkers in two spatial dimensions. To mimic turnover of 
axons, whole paths are removed and new walkers are injected with specified rates. In the simplest 
version of the model, we use strongly adhesive short-range inter-axon interactions that are identical 
for all pairs of axons. We generalize the model to adhesive interactions of finite strengths and to 
multiple types of axons with type-specific interactions. The dynamic steady state is characterized 
by the position-dependent distribution of fascicle size and fascicle composition. With distance in the 
direction of axon growth, the mean fascicle size and emergent time scales grow monotonically, while 
the degree of sorting of fascicles by axon type has a maximum at a finite distance. To understand 
the emergence of slow time scales, we develop an analytical framework to analyze the interaction 
between neighboring fascicles. 

PACS numbers: 87.19.1p, 05.40.-a, 05.40.Fb, 87.19.1x 



I. INTRODUCTION 

Reaction-diffusion phenomena arise in diverse fields 
such as physical chemistry jlj or developmental biol- 
ogy 0. In certain reaction-diffusion systems, the process 
of path aggregation occurs, in which preferred paths of 
the diffusing elements are established and evolve in time. 
The path aggregation process is found in diverse realm of 
nature, e.g., in formation of insect pheromone trails 
and human walking trails (g, |jj , in aggregation of trails 
of liquid droplets moving down a window pane, in river 
basin formation [1, Q, etc. 

One class of the mathematical models in which path 
aggregation processes have been studied is the active- 
walker models 0, fl0| in which each walker while moving 
through the system changes the surrounding environment 
locally, which in turn influences the later walkers. The 
ant trail formation is an example of such a process 0, 0] ■ 
An ant leaves a chemical trail of pheromones on its path 
which the other ants can sense and follow. Evaporation 
of pheromone leads to an aging of these trails. Similarly, 
the mechanism of human and animal trail formation is 
mediated by the deformation of vegetation that generates 
an interaction between earlier and later walkers @, Q. 
This deformation, and therefore its impact, decays con- 
tinuously with time la- 
in a recent Letter we analyzed the dynamics of 
path aggregation using a simple model that belongs to 



* Current address: FOM Institute AMOLF, Science Park 104, 
1098 XG, Amsterdam, The Netherlands; Electronic address: 
d. chaudhuri® amolf. nl 

^Electronic address: peterphysikOgmail.com 
t Electronic address: zapotocky@biomed.cas.cz 



the class of active walker systems discussed above. In 
contrast to the active-walker models, in our model the 
individual paths do not age gradually, but rather main- 
tain their full identity until they are abruptly removed 
from the system. This particular rule for path aging was 
chosen to allow application of our model to the process 
of axon fasciculation (formation of axon bundles [13). 
which we discuss below and more in detail in Sec. ITT] 

In order to develop neuronal connections, sensory neu- 
rons born in peripheral tissues project their axons (long 
tubular part of the neuron cell that conducts electrical 
excitations) towards target regions in the brain. Fre- 
quently, multiple axons come together to form axon fas- 
cicles, and may sort according to the cell type of the neu- 
ron to which the axon belongs. This fascicle formation 
and sorting can be driven by inter-axon interactions lead- 
ing to, e.g., a pre-target spatial map in the mammalian 
olfactory system [pi - fla ]. 

In our model, the axons are represented as directed 
random walks in two spatial dimensions. In Ref. 
we formulated and analyzed the simplest version of the 
model, in which all axons belong to the same type and 
have strong adhesive interactions, so that each newly 
growing axon encountering an existing fascicle will join 
the fascicle and never detach. In the presence of axon 
turnover (aging of the paths), a steady state character- 
ized by a distribution of fascicle sizes is eventually 
established. The focus of Ref. [ll[ was on the analysis of 
the surprisingly long time scales that emerge from this 
simple dynamics. 

In the current paper, we significantly extend this theo- 
retical analysis. We develop an analytical description of 
the dynamics of two neighboring fascicles, and show how 
the slowest mode of their interaction gives raise to the 
slow time scales observed in Ref. [ll[ . We also systemat- 



ically discuss the limited analogies that can be made be- 
tween our 2-dimensional model and 1-dimensional models 
of particle coalescence [Itij . aggregation [13], and chip- 
ping [lj| [l9|] . These analogies are useful for the under- 
standing of stationary quantities of our model such as 
the distribution of fasicle sizes and the distribution of 
inter-fascicle separations. 

The main contribution of the current paper, however, 
is to generalize the previous model of Ref. [1JJ to attrac- 
tive interactions of finite strength (so that detachment 
of axons from fascicles is possible) and to multiple axon 
types with type-specific interactions. Such a generaliza- 
tion is necessary to allow the biological application of the 
model. 

In the following section we give a detailed biological 
motivation for the model we consider. In Sec. IIIII we 
introduce the model and the Monte Carlo (MC) simu- 
lation scheme that we use to investigate its properties 
numerically. Followed by this, in Sec. IIV1 we give a brief 
overview of guiding concepts that will recur in the rest of 
the paper. In Sec. El we present a detailed analysis of the 
system containing axons of a single type that follows the 
"always attach, never detach" rule. We extend the nu- 
merical results of Ref. [HI for the properties of the steady 
state and for the emerging time scales. We review the 
analytical framework of single-fascicle dynamics, devel- 
oped in Ref. fTlj j . and significantly extend it by deriving 
results for the interaction dynamics of two neighboring 
fascicles. In Sec. I VII we numerically study the effects of 
non-vanishing detachment rates of axons from a fasci- 
cle. In Sec. IVIII we discuss some limited analogies of our 
model to one dimensional aggregation and coalescence 
processes. In Sec. I Villi we discuss the sorting of fascicles 
by axon types in a system containing two types of axons, 
the simplest manifestation of a mixed population of ax- 
ons. Finally, we provide a summary of main results in 
Sec. HXI and conclude in Sec. [X] by discussing the outlook 
for biological applications of our model. 



II. BIOLOGICAL MOTIVATION 

Sensory neurons located in peripheral tissues connect 
to more central locations of the nervous system via, ax- 
ons [2pj . During the development of an organism, axons 
of newly maturing sensory neurons must establish con- 
nections to the proper location. Axon growth is initiated 
at the soma (main cell body) of each neuron, and pro- 
ceeds with a typical rate of extension 100/xm/h [2lJ. The 
direction of growth is controlled by the dynamic growth 
cone structure at the tip of the axon. The growth cone 
probes the environment in its vicinity, and can detect 
gradients of spatially distributed chemical signals. In 
the absence of strong directional signals, the path of the 
growth cone is highly stochastic [22|, [23|, while in the 
presence of appropriate guidance cues, the direction of 
motion becomes strongly biased. The overall direction of 
axon growth may be guided by spatial gradients of chem- 



ical cues generated by the target. A number of distinct 
molecular guidance cues that influence neuronal devel- 
opment have been identified in recent years 0, , and 
the response of the growth cone to graded cues has been 
studied theoretically [26- 28] . In this work, we do not di- 
rectly model the axon guidance by graded chemical cues, 
but subsume their influence into the setup of our model 
by giving all axons a common preferred growth direction. 

In this article, we study the collective effects that arise 
from direct local interactions among the growing axons. 
When such interactions are attractive, the growth cone 
of a newly growing axon tends to follow the tracks (i.e., 
the axon shafts) of older axons. The strength of this in- 
teraction is governed by the type and expression level of 
the relevant cell adhesion molecules [U [2{| . The result- 
ing dynamics can lead to selective formation of fascicles 
of axons [30U32T ] , a common and essential phenomenon in 
the developing nervous systems. 

An additional important aspect included in our model 
is that of neuronal turnover. During development, a sig- 
nificant portion of sensory neurons with fully grown ax- 
ons may die, and be replaced by younger sensory neu- 
rons which attempt a new connection to the brain. For 
example, up to 80% of retinal ganglion cell axons are 
lost dur ing, the development of the visual system in the 
cat |33l . l34j . In the mammalian olfactory system, both 
neuronal birth and death persist throughout the life of 
the animal, leading to a dynamical steady state pattern 
of connectivity. In particular, the average lifetime of an 
olfactory sensory neuron in the mouse is of the order of 1- 
2 months [Hj], which is less than one tenth of the mouse's 
lifespan. 

To motivate the introduction of multiple types of ax- 
ons into our model, we now briefly discuss the intricate 
connectivity pattern of the mammalian olfactory system, 
which implements the sense of smell. In the mouse, the 
adult nasal epithelium contains approximately 10 6 olfac- 
tory sensory neurons, which send their axons through the 
olfactory tract to the olfactory bulb in the forebrain. Re- 
markably, the sensory neurons belong to approximately 
1200 distinct types 0, [36j], and the axons of each type 
connect to a distinct neuronal structure, a glomerulus, 
on each olfactory bulb [HI, [37} ■ Such precise connectiv- 
ity is fully established only after several turnover periods, 
while in newborn mice, split glomeruli and glomeruli that 
mix several axon types are often observed [35l , l38j . 

In olfactory sensory neurons, elegant genetic analysis 
shows that the axonal type is determined by the ex- 
pression of a specific odorant receptor gene [H, l39M4l| . 
Physiological experiments on mice show that the expres- 
sion of specific types of cell adhesion molecules, that dic- 
tates the strength of adhesive forces between axons, is 
strongly correlated with this axonal type [40j, |42j . A wide 
range of strengths of interactions between axons may be 
generated through combinatorial expression of multiple 
types of cell adhesion molecules. 

In Fig. Q] we show configurations of olfactory axons as 
observed in in vivo filj ] and in vitro [43j experiments. 
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FIG. 1: (Color online) (a) Axons of olfactory sensory neurons of a specific type (M71) growing in the surface layer of the mouse 
olfactory bulbs. Scale bar = 500 /im. The axons emerge (top) from behind the olfactory bulbs, and grow towards the bottom. 
Note fasciculation as well events of detachments of axons from fascicles. Figure adapted from Ref. [43]. (b) Axon growth (top 
to bottom) and fasciculation observed in explant culture of rat olfactory epithelium on a laminin-coated coverglass. Scale bar 
= 100 /im. Axon type is not distinguished. Figure adapted from Ref. [431 ]. 



Fig. [TJa) shows axons growing in the surface layer of 
the left and right olfactory bulbs of a genetically modi- 
fied mouse (Fig.l(L) of Ref. [4l|). Only axons belonging 
to one type of olfactory sensory neurons (expressing the 
M71 receptor gene) are labeled; axons of other types are 
present but not visible. The axons progressively fascicu- 
late and terminate in a glomerulus visible in the center 
of each half-image. Fig.QJb) shows fasciculation of axons 
growing from an explant of the rat olfactory epithelium 
(Fig. 7(a) of Ref. [43]]). In this case, the fluorescent la- 
beling does not distinguish the axonal type, and (with 
a high probability) the visible axons belong to multiple 
types. 

Our model aims to provide a quantitative framework 
for evaluating the contribution of axon-axon interactions 
to the formation of patterns described above. The pres- 
ence of turnover and multiple axon types in our model 
distinguishes our work from previous theoretical studies 
of axon fasciculation [3, EEJ . Our implementation of the 
individual axon dynamics is particularly simple, to al- 
low us to concentrate on collective effects arising from 
interactions within a population of axons. 



III. MODEL AND NUMERICAL 
IMPLEMENTATION 

A. Setup and interactions 

In our model, each growing axon is represented by a di- 
rected random walk in two spatial dimensions (Fig.[2ja)). 
The random walkers (representing the growth cones) are 
initiated at the epithelium (y = 0, random even x) with 
a birth rate a, and move towards the bulb (large y) with 
constant velocity v y = 1. In the case of multi-type sys- 
tems, a type is assigned to each newly initiated random 
walker (specifically in the simulations of Sec. IVIII1 the 



type is decided randomly with equal probability for each 
of the two types). The trail generated by a random 
walker (growth cone) is regarded as an axon shaft. A 
forward moving directed random walker (growth cone) 
interacts with trails (axon shafts) of other walkers. In 
the numerical implementation on a tilted square lattice, 
at each time step the growth cone at (x, y) can move to 
(x - l,y + 1) (left) or (x + l,y + 1) (right). The proba- 
bility P{l,r\ to move left/right is evaluated based on the 
axon occupancies at the (x±l, y + 1) sites and their near- 
est neighbors (x ± 3, y + 1) (see Fig. [21a))- At a given y, 
two axons are considered to be part of the same fascicle if 
they are not separated by any unoccupied sites (i.e., they 
are not separated by more than two lattice spacings). 

We assume a short-range attractive interaction be- 
tween each growth cone and the close-by axon shafts. 
The range of interaction (two lattice spacings in our 
model) corresponds biologically to the range of extension 
of sensory filopodia from the growth cone (of the order 
of 10/im). The attractive interaction is mediated by cell 
adhesion between the growth cone and the axon shafts. 
We assume that the interactions are additive and type- 
specific. For a given growth cone, the model assumes a 
weak nearest neighbor attraction E Q < if the neigh- 
boring axon shaft is of a different type and a stronger 
attraction Eh{< E Q ) if the neighboring axon shaft is of 
the same type. In each time step, a growth cone at (x, y) 
attempts a Monte Carlo move to the left (x — l,y + 1) 
and to the right (x+l,y + l) with probabilities 1/2. The 
moves are accepted with probabilities 

Pl = min[l,exp(— 8 Ei)] 

(or pr = min[l, exp(— 5E r )]) where 

SEi = [rih(x-3,y + l)-rih(x + 3,y + l)]E h 
+ [n (x — 3,y + 1) — n (x + 3, y + 1)]E 
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FIG. 2: (Color online) (a) Interacting directed random walks 
on a tilted square lattice. A random walker (+) represents 
a growth cone. For one walker, the possible future sites (□) 
and their nearest neighbors (o) are marked. The trail of a 
walker (line) models an axon shaft, (6) A typical late-time 
configuration (t — 25T) of a system of axons belonging to 
two different types, r (red) and b (blue). The strength of the 
homotypic interaction is Eh = —4 and that of the heterotypic 
interaction is E = —0.1. The mean numbers of r and b 
axons at y = are Nq = Nq = 50 and the system size is 
L = 800. (c) A typical late-time configuration (t = 25T) 
in a system with a single type of axons undergoing energy- 
minimizing dynamics with system parameters L — 800, Ao = 
100. For the fascicle identified at y = 600 (arrow), D indicates 
its basin and E is the inter-basin free space. 

(SE r = —8 Ei) and the occupancy number rih denotes 
the number of axons belonging to the same type as the 
growth cone, while n Q is the number of axons of other 
types. Notice that in calculating the difference in energy, 
the occupancy of the positions (x ± l,y + 1) does not 
appear, as their contributions to the energy cost mutually 



cancel. Periodic boundary conditions are used in the in- 
direction. In this kinetic MC scheme, we used parallel 
updates of all random walkers in each MC step. 

We now clarify the relation of the general model de- 
scribed above to the model studied by us in Ref. I n 
this simple version of the model, all axons belonged to 
the same type and the interaction between them was gov- 
erned by the "always attach, never detach" rule, which is 
the "zero temperature" version of the above-mentioned 
general MC scheme, i.e., the dynamics was a pure energy 
minimization process: Pl = I (pu = 1) when SEi < 
(SE r < 0); pl = Pr = 1/2 in all other cases. This is 
in contrast to the "finite temperature" dynamics of the 
general model, in which there is a non-zero rate for the 
detachment of growth cones from fascicles. Ref. [ll[ used 
sequential updates, in contrast to parallel updates used 
in this paper. An additional difference is that in Ref. [TTI j . 
the interaction was not assumed to be additive, i.e., the 
strength of interaction of a growth cone with a fascicle 
did not depend on the number of axons in the fascicle. 

Note that two fasciculated axons run parallel to each 
other with their separation in x-direction restricted 
within two lattice spacings, the interaction range. Thus 
the typical width of a fascicle containing multiple axons 
remains 2-3 lattice spacings in our model. We do not 
implement any on-site repulsion, i.e., axons are free to 
grow on top of each other. 

In our model, we do not consider any relaxation dy- 
namics of axon shafts. This corresponds to the assump- 
tion of strong adhesion of axons to the substrate, so that 
the line tension on axons can not straighten out the local 
curvatures. 



B. Turnover 

To capture the effect of neuronal turnover, each 
random walker is assigned a finite lifetime 9 from 
an exponential distribution of lifetimes H(9)d0 — 
±exp(-6/T)d6 with mean (0) = J™ 9II{8)d8 = T. 
When the lifetime expires, the random walker and its 
entire trail (i.e., the whole axon) is removed from the sys- 
tem. The mean number of axons in the system reaches 
the steady-state value N(y) = Noexp(—(3y), where (3 = 
1/T is the mean death rate per axon and the steady- 
state occupancy at y — is Ao = a/f3. In the simula- 
tions, we use T = 10 5 time steps, and restrict ourselves 
to y < T/10. The birth rate a is chosen so as to obtain 
the desired number of axons No, or equivalently, the de- 
sired axon density p = Nq/L (ft = 1/2 implies an average 
occupancy of one axon per site), where L is the system 
size in cc-direction. As we will show in detail, the time 
scales needed to achieve the steady state of fascicle size 
distribution can be very long compared to the time scale 
T needed to achieve the steady state value of the total 
number of axons. A typical late-time configuration for a 
system of axons involving two distinct axonal types with 
type-specific interactions is shown in PigJ2j6). 
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C. Parameters 



TABLE I: Parameters of the model 



In this section we briefly discuss the biological meaning 
and physical values of the parameters in our simulations. 
In our model the interaction range in x- direction is cho- 
sen to be 2 lattice units. Since we assume only contact 
interactions, this range of 2 lattice units corresponds to 
the length of a nlopodium which typically is 10 /im. Thus 
the lattice spacing in x-direction Ax = 5 /im. 

The time step At in our model needs to be large enough 
to allow the growth cone to integrate a signal and react 
to it. Ref. [23] suggests this time scale to be of the order 
of tens of seconds; we choose At = 60s. This corresponds 
to a diffusion constant (in x- direction) (Ax) 2 /2At = 
12.5 /xm 2 /minute which compares well with in vitro ob- 
servation for short time scales up to tens of minutes [22j ■ 
Note that with this choice of At the mean lifetime of 
an axon T = 10 5 At used in our simulations corresponds 
to 69.4 days, quite typical of axons of mouse olfactory 
sensory neurons [35| . 

The lattice unit in y-direction is now chosen to give 
a reasonable growth velocity v y . Choosing Ay = 1 /zm, 
we have v y = 60/im/hour, which is a typical value for 
growing axons of sensory neurons [2l| . 

Note that with the above mentioned choices, Fig. &b) 
corresponds to a system size of 800 Ax = 4 mm in IE- 
direction and 10 4 Ay — 10 mm in the y- direction. These 
dimensions are comparable to the size of the olfactory 
bulb in mice [36j |. 

The effective interaction energies Eh and E a should 
be chosen to match the observed rates with which ax- 
ons detach from fascicles. We introduce the quantity 
TTd = exp(E)/Ay that expresses the rate of detachment 
of one axon per unit length of a two-axon fascicle given 
that the two axons interact via E — Eh or E a . Thus a 
growth cone interacting with a fascicle of n axons will fol- 
low the fascicle over the mean distance L y — n~ d n before 
it detaches. It is not straigthforward to use published 
experimental images to deduce E as usually, the size of 
the fascicle is not known, and the location at which a 
growth cone first attached to the fascicle is not recorded. 
The observed typical distance L y varies widely depend- 
ing on the specific neural system, ranging from tens of 
/j,m to centimeters. In our simulations, we use the range 
of homotypic interaction strength Eh = —4 to -1, which 
corresponds to detachment rates TTd = 0.02 to 0.37 /im -1 . 
In Table U we list the meaning and values of the param- 
eters used in our model. 



IV. OVERVIEW 

We show in this paper that simple directed growth of 
axons that interact via a short-range attraction leads to 
reliable formation of axon fascicles, in absence of any ex- 
ternal chemical guidance cue. Once a fascicle is formed 
its position does not move appreciably, however, the fas- 
cicle size (number of axons present in the fascicle) fluctu- 



Symbol 


Meaning 


Value 
(simulation) 


Value 
(physical) 


Ax 


Lattice spacing 
in ir-direction 


1 


5 pm 


Ay 


Lattice spacing 
in y-direction 


1 


1 pm 


At 


Time step 


1 


60s 


T 


Mean axonal lifetime 


10 s 


69.4 days 


N 


Mean number of 
axons at y = 


50 to 200 


50 to 200 


L 


System size 
in 3>direction 


100 to 800 


0.5 to 4 mm 


E h 


Homotypic 
interaction strength 


-4 to -1 


detachment rate 
TTdh = 0.02 pm' 1 
to 0.37 pm' 1 


Eo 


Heterotypic 
interaction strength 


-0.1 


detachment rate 
n do = 0.9 fim -1 



ates. The turnover of individual axons generates a slow 
dynamics of reorganization of fascicles at a fixed y-level. 

In the simplest case (Sec.fV]), the system contains only 
a single type of axons which grow and form fascicles via 
an energy-minimizing dynamics (strong inter-axon inter- 
action), so that once attached the axons do not leave a 
fascicle. For this case one can uniquely assign a basin of 
each fascicle at any specified y-level (see Fig. [H[c)). The 
basin size D of a fascicle is the interval at the level y = 
between the right-most and left-most axons belonging to 
the fascicle (Fig. Hfc)). Any axon growing from within 
this basin must contribute to the fascicle size unless it 
dies before reaching the specified y-level. Thus the aver- 
age number of axons that survives in a fascicle at level 
y (in the steady state) is n = D pexp(—f3y). The axons 
initiated at the opposite edges of the basin are expected 
to meet each other in y ~ (D/2) 2 steps of random walk 
in x-direction. Therefore, one obtains the mean-field pre- 
diction for the mean fascicle size n(y) ~ 2y x l 2 p exp(—/3y) 
up to y ~ (i/2) 2 , where complete fasciculation occurs, 
i.e., n = N(y). 

In a system with finite detachment rates (Sec. IVI[) . 
however, growing axons can leave one fascicle and at- 
tach to another. Thus the fascicle basins overlap and 
axons introduced in the basin of one fascicle can end up 
in a different fascicle. Still the mean-field estimate of 
the increase of mean fascicle size with y, shown above, 
turns out to remain approximately valid (Sec. IVI"S]) . For 
higher detachment rates (weaker interactions), the pref- 
actor of the power-law growth is reduced, corresponding 
to smaller fascicles. In the limit of extremely weak in- 
teractions, each axon would grow independently of the 
others, and no fasciculation is possible. This overall pic- 
ture remains intact even for systems having multiple axon 
types. The dynamic steady state is characterized by a po- 
sition (y) dependent distribution of fascicle sizes which 
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shows a scaling law. The peak of the distribution shifts 
towards larger fascicle sizes at higher y-levels (Sec. IV Bl 
and Sec. EE1]). 

In a steady state configuration at fixed time t (such 
as in Fig. HJc)) the axon fasciculation with increas- 
ing y may be formally viewed as the evolution of 
an one-dimensional diffusion-aggregation or diffusion- 



coalescence process [46|. These limited analogies can 



be used to approximately understand steady state prop- 
erties like distribution of fascicle sizes and inter-fascicle 
separation (Sec. IVIljl . However, the full dynamics of our 
model has no simple one-dimensional counter-part [ill ]. 
A projection of the dynamics onto one dimension would 
involve complicated long-time correlations between the 
random walkers. 

The dynamics of fascicle reorganization can be charac- 
terized by the slow approach to steady state, or by the 
steady state auto-correlation time for the mean fascicle 
size. These time scales grow with y and may reach values 
orders of magnitude larger than T. In absence of detach- 
ment, the slowest mode of fascicle reorganization occurs 
via partial exchange of neighboring basins (Sec. IV D|) . 
In the presence of detachment, in addition, the basin 
of one fascicle can easily drain to another. Thus the 
time scales decrease with decreasing inter-axon attrac- 
tion fScc lVTAl 

In systems containing multiple types of axons with 
type-specific interactions, we evaluate the degree of sort- 
ing S that quantifies the type-wise purity of the local 
environment of axons (Sec. I VIII j) . At steady state, S(y) 
shows a non-monotonic variation, with a maximum at 
some intermediate y. The position and the value of the 
maximum depend on system parameters like the mean 
density of axons p and the interaction strengths. This 
non-monotonicity is due to the attractive heterotypic in- 
teraction which merges mid-sized, relatively pure fasci- 
cles to form large impure fascicles. 



SINGLE TYPE OF AXONS, NO 
DETACHMENT 



In this section we analyze the collective behavior of 
axons belonging to a single type following the energy- 
minimizing "always attach, never detach" rule. This 
model has been investigated in detail in a previous publi- 
cation In this section, we extend the numerical and 
analytical results of Ref. [11]. In the simulations, we use 
a modified implementation of the Monte Carlo update 
rules. In contrast to Ref. [llj , the strength of interaction 
between a growth cone and a fascicle is assumed to be 
proportional to the number of axons present in the fas- 
cicle. Another difference is that we use parallel updates, 
instead of sequential updates used in Ref. ll| . We there- 
fore include a comparison to the main results we reported 
in Ref. [ill, to show that these are not altered. 




FIG. 3: (Color online) Approach to the steady state in the 
L = 400, No = 200 system at representative j/-levels indicated 
in the legend. The mean fascicle size (n(t;y)) (averaged over 
10 3 realizations ) approaches n oa {y) as t — > oo. The data set 
labeled as p — 1/8 is from the L = 400, No = 50 system, 
collected at y — 10 4 . Inset: Fitting of rtoa — (n(t,y)) to a 
function f(t) — pexp(—/3t) + q exp(— t/r ap ) shown in a semi- 
log plot. The data is the same as shown in the main figure 
at y = 10 3 . The fitting parameters are rioo = 42.94 ± 0.07, 
p = 13.98 ± 0.24, q = 7.43 ± 0.05, and approach-to-steady- 
state time scale r ap = (294 ± 6)T. We used the Marquardt- 
Levenberg algorithm for nonlinear least-squares fitting as im- 
plemented in gnuplot version 4.4. 



Approach to steady state: mean fascicle size 
and time scale 



A typical late-time configuration for a system with 
L = 800 and N = 100 (density p = N /L = 1/8 at 
y = 0) is shown in Fig. H][c). With increasing y, the ax- 
ons aggregate into a decreasing number of fascicles. The 
number of axons in a fascicle is referred to as the fasci- 
cle size n. At steady state, the mean fascicle size n at 
level y may be estimated as n ~ Ipy 1 / 2 exp(— fly) up to 
y ~ (L/2) 2 , where complete fasciculation n = N(y) is 
expected [llj . 

The measured mean fascicle size, obtained by averag- 
ing over all the existing fascicles at a given y (Fig. [3]), 
grows with time asn — rioo — pexp(— flt) — gexp(— t/r ap ), 
where T ap (y) defines the time scale of approach to the 
steady state value n^y). The same behavior was ob- 
served earlier in simulations reported in Ref. [11]. The 
semi- log plot in the inset of Fig. [3] shows clearly the slow 
exponential approach to the steady state mean fascicle 
size. Using the above-mentioned double-exponential fit- 
ting we extract the time scale r ap and the steady-state 
mean fascicle size rioa at all the y-levels. 

The approach-to-steady-state time scale r ap increases 
with y. T ap can exceed the mean axon lifetime T by orders 
of magnitude (Fig. [3J and QJ. Note that T ap is longer in 



11] 



a system with larger density of axons p (Fig. @|. Ref 
discussed this point in detail. Further, asymptotically 
in y, we find n^, — c + 2py b exp(— fly), with b w 1/2 
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FIG. 4: (Color online) Approach-to-steady-state time scale 
Tap as a function of y at different axon densities. Tap is ex- 
tracted by a fitting procedure as described for the inset of 
Fig. [3] The time series are collected over t = 500 T and 
averaged over 10 3 realizations . All the data were collected 
for a system of No = 50 axons, with varied system sizes, 
(i) L = 100 = 1/2), (ii) L = 150 (p = 1/3), and (iii) 
L — 200 (p = 1/4). Correlation time t c : The last data set 
shows the steady state correlation time r c (in units of T) 
for a L = 100, No = 50 system calculated from the corre- 
lation function c(t) = (n(t)n(O)). c(t) is evaluated by us- 
ing the time series of n(t) collected between t = 200T and 
2 x 10 4 T, and averaged over 30 realizations . A fitting of 
c(t) = p + qexp(— fit) + r exp(— t/r c ) allows us to extract r c at 
different y-levels. Fitting errors in T ap and r c are within 5% 
(see the caption to Fig. [3] (inset) for the fitting procedure and 
estimate of error in T ap ). The thick solid line shows a power 
law y 2b with 6=1/2. 



(Fig. [5]) - in good agreement with the mean-field predic- 
tion (Sec. TQ- 

Impact of interaction range: To test the impact of the 
range of inter-axon interaction, we have simulated a sim- 
ilar system with purely "contact" interaction, i.e., the 
interaction range is taken to be zero. With this reduc- 
tion in the range of interaction, we find that the emerg- 
ing time scales decrease. For instance, for a system of 
L = 800 and N Q = 100 the approach-to-steady-state time 
becomes r ap < 10T. Thus an increase in the range of in- 
teraction increases the emerging time scales. The steady 
state distribution of fascicle sizes shows the same scaling 
behavior as in the case of nearest-neighbor interaction 
discussed in the following. 



B. Steady state 

The steady state is characterized by the stationary dis- 
tribution of fascicle sizes P s (n, y), defined as the number 
of fascicles of size n at level y. 




100 1000 10000 

y 



FIG. 5: (Color online) Time-asymptotic fascicle size rioo 
for systems with No = 50 axons and system sizes L = 
100, 150, 200, 400 (corresponding to p = 1/2, 1/3, 1/4, 1/8 
respectively) as a function of y. The time series are collected 
over t = 500 T and averaged over 10 3 realizations . The sub- 
tracted mean fascicle size n s = (rioo — c) exp(/3y)/2p corrected 
for the finite axon- lifetime T = 1//3 is plotted as a function of 
y. The offset fascicle size c is treated as a fitting parameter, 
c = 8.17, 3.48, 1.56, 0.11 for L = 100, 150, 200, 400 respec- 
tively. Fitting errors in n B are within 3%. Data collected 
at the various densities collapse onto a power law y b with 
b = 1/2. The largest system shows the widest power law 
regime. 

1. Steady state: scaling regime 

For a system with L = 800 and -/V = 100, P s (n,y) is 
shown at a series of y-levels in Fig. [5] Within the range 
y = 10 3 — 10 4 all data collapse onto a single curve after 
appropriate rescaling (Fig. [5]). This data collapse implies 
the scaling law [llj 

P s (n,y) = (n(y))- r cb(n/(n( y ))) (1) 

with r = 2.1 and the scaling function 4>{u) — 
J\fuexp{—vu — \u 2 ). Note that the steady state aver- 
aged fascicle size (n(y)} is a quantity equivalent to the 
asymptotic n^y) discussed in the previous subsection. 

The scaling law in Eq. [T]can be justified starting from 
the assumption of homogeneity of fascicle size distribu- 
tion P s {n,Xy) = \~ p P s {\~ q n,y), with the exponents p 
and q undetermined at this stage. Noting that the mean 
number of axons N(y) = J dnn P s (n,y) and the mean 
number of fascicles B(y) = J dnP s (n,y), the homogene- 
ity condition leads to the relations N(Xy) = \~' p+2q N(y), 
B(\y) = \~ p+q B(y). Since by definition the mean fas- 
cicle size (n(y)) = N{y)/B{y), (n(\y)} = \ q (n(y)). In- 
voking the mean-field prediction (n(Xy)) — X b (n(y)) with 
6=1/2 (Sec. E] and El A), we find q = b. If N(y) were 
independent of y, we would have had p = 26. However, 
in fact N(y) = Noexp(—/3y). In the region j3y < 1, 
we can write p — 2b + S with <5 w /3y/]ny <C 26. 
Note that the relation P s (n,Xy) = \~ p P s (\~ b n,y) can 
be recast in the form P s (n,y) — (n)~ r (j)(n/(n)) where 
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sharply peaked at n = N(y). To demonstrate this fact 
we take a system of small size and high density p = 1/2 
(L = 100 and No = 50). Within the scaling regime 
(y 500) the distribution function maintains the scal- 
ing form uexp(— vu — Xu 2 ). However, at higher y-levels 
the distribution becomes bimodal with a new maximum 
appearing, characteristic of the complete fasciculation. 
This shows a coexistence of two preferred fascicle sizes. 
Finally, at y ~ (L/2) 2 the whole weight of the distri- 
bution shifts to this new maximum and the distribution 
becomes unimodal again (see Fig. [7]). In systems with 
larger L, the regime of coexistence shifts towards higher 
y- levels. 



FIG. 6: (Color online) Steady state distribution of fascicle 
sizes P s (n, y) (averaged over 10 4 realizations and the time 
interval 10T < t < 25T) for the No = 100, L = 800 sys- 
tem at y- levels indicated in the legend. Inset: A scaling 
with B = l/(n) and A = (n) 2 ' 1 collapses all data obtained 
for y = 1585, 1995, 3162, 5012, 6310, 7943, 10 4 onto a sin- 
gle curve 4>{u) = Afu exp(— vu — Am 2 ) with u — n/(n) and 
N = 274, v = 0.78, A = 0.45. 



y = 1000 




1 15 30 45 60 



n 

FIG. 7: (Color online) Fascicle size distribution at large y 
for a system of L — 100 and iVo = 50. The data were av- 
eraged over 400T < t < 500T and 10 3 realizations . The 
single-peaked distribution characteristic of the scaling regime 
(y = 500) crosses over to a distribution peaked near com- 
plete fasciculation n = N(y) = 48.77 at y = 2500 [= (L/2) 2 ] 
through a coexistence regime showing a double maximum (at 
y = 1000). 

r = p/b = 2 + 5/b > 2, in agreement with Eq. Q] As 5 is 
y-dependent, the scaling of P s (n,y) is only approximate. 

2. Steady state: crossover to complete fasciculation 

The steady state distribution changes its shape dras- 
tically beyond the scaling regime. Near y = (L/2) 2 , 
on an average all the axons are expected to collapse 
onto a single fascicle, thereby generating a distribution 



3. Steady state: correlation time 

The dynamics in the steady state is characterized by 
the auto-correlation function for the mean fascicle size 
n(t) at a fixed y-level: c(t) = (h(t)h(0)) which fits to 
the form p + qexp(—(3t) + r exp(— t/r c ) (as in Ref . (TIT] ) . 
The correlation time t c increases with y and significantly 
exceeds the axon lifetime T. We show this behavior for 
systems at p = 1/2 (L = 100, iV = 50) in Fig. [U This 
shows a regime of approximate power law growth of the 
time scale t c ~ y 2b with b m 1/2. 



C. Effective single- fascicle dynamics at fixed y 

In this subsection, we review our analytical results 
from Ref. [11]. The following subsection presents new 
results for the time scales arising from the interaction of 
two neighboring fascicles. 

The concept of effective single-fascicle dynamics has 
been introduced in Ref. [ill] . The dynamics of a mean 
fascicle at level y with n(t) axons can be viewed as a 
stochastic process with gain rates u+(n) (for transitions 
n — > n.+ 1) and loss rates u~ (n) (for transitions n — > n — 
1). A fascicle loses axons only by the death of individuals, 
therefore, U-(n) = fin fill ]. 

In absence of detachment events, any axon introduced 
within the basin (size D) of a fascicle (see Fig.[2][c)) can 
not escape the fascicle. Moreover, some of the axons born 
in the neighboring inter-basin gaps (size E) eventually 
join the fascicle under consideration. These two processes 
contribute to u + (n). As was shown in Ref. [ll|, the time 
series of D(t) and n(t) tend to co-vary. Thus treating 
the dynamics of D as slave to n, we get a form u + = a + 
bn [ll|. Note that the basin size D can not exceed 2y or 
L, and D > 2y x / 2 occurs with low probability. Therefore 
a saturation of u+(n) is expected for large values of n. 
The measured average gain and loss rates u± (n) obtained 
from our current simulations agree with the functional 
forms u + (n) = a + + b + n — c + n 2 and u_(n) — j3n (data 
not shown). The quadratic correction to linear growth 
captures the saturation of u+(n) at large n. 
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The master equation of the growth-decay process for 
the effective single fascicle of size n at level y may be 
written as 



u + {n- l)P{n- l,t)+u-(n + l)P(n + l,t) 
[u+{n)+U-(n)]P(n,t), (2) 



P(n,t) 

for n > 1. For the boundary state (n = 1 



P(l,t) = J+(y) + u.(2)P(2,t) - [u+(l)+«_(l)]P(l,t) 

where J+{y) represents the rate with which new single 
axons appear between existing fascicles at y. 

The solution of the master equation at steady state 
was derived in Ref. [ll[ and has the form, 

pP s (n,y) = J+(y) rf exp[-i(n - 1) - «(n- l) 2 ], (3) 

where 7 = a + /(3 — 1 and I = f — b + / (3 and k = c+/2{3. 

From the master equation one can estimate the 
approach-to-steady-state time r ap and the correlation 
time at steady state r c (llj . The correlation time t c 
for the fascicle size n, near the macroscopic stationary 
point n s [u + (n s ) = u_(n s )] can be expressed [47| as 
t c = l/(«'_(n a ) - <(n s )) = 1/08 - 6+ + 2c+n s ). Un- 
der the linear approximation of u + (n) = a + + 6 + n the 
approach-to-steady-state time scale for the average fasci- 
cle size (n) can be written as r ap — l/(/3 — b+) [43]. 

Further, the mean lifetime of fascicles can be defined 
as Tf = [f^° P s (n,y)dn]/J + (y) and is evaluated to ob- 
tain [11 



(T/2k) 1 - ^V^e^ (I - 2re)erfc(£/2 N /K)) /2>/k 



Notice that the above derivations of the time scales 
already involved the numerical observation u + (n) — 
a + + b + n — c + n 2 . Using the y-dependence of b + and 
c+ obtained from numerical simulations, we found power 
law growth of the time scales, t c ~ y b , r ap ~ y b and 
T f ~ y 2b with 6 fa 1/2 [ill]. 

In the following subsection, using a purely analytical, 
deterministic treatment of the dynamics of two neigh- 
boring fascicles, we show that a time scale growing as 
y 2b emerges due to an exchange of basin size between the 
fascicles. 



D. Effective dynamics of two interacting fascicles 

In this section we analytically explore the effective 
dynamics of two neighboring fascicles. Under the "al- 
ways attach, never detach" rule the basins of neighbor- 
ing fascicles do not overlap. Three dynamical variables 
characterize the dynamics of the fascicles: (i) the num- 
ber of axons ni present in the fascicle, (ii) the basin 
size Di and (iii) the separation Ei between the i-th and 
(i + l)-th basin (Fig. (He)). We explore the dynamics 
at a fixed y. We consider the thermodynamic limit of 
large Nq and L (with a fixed density p = Nq/L), and 




FIG. 8: (Color online) Illustration of basin size exchange, i.e., 
the slowest mode in the effective dynamics of two fascicles 
(Sec. IV D)) . Two neighboring fascicles i and i + 1 are shown 
along with their corresponding basins Di and -Di+i, and inter- 
basin gaps Ei and -Ei+i. The exchange of the boundary axon 
(dashed green line) between the two fascicles corresponds to a 
relaxation mode with time scale r ~ y 2b (see the main text). 
This increases the basin size Di+i at the cost of Di, leaving 
the gap size Ei unaltered. 



express all the length scales in units of L and number 
of axons in units of Nq. Thus the reduced variables are 
Vi = n i/No, \i = Di/L, ei = Ei/L, and the reciprocal 
system size a = 1/L has the meaning of a lower cut-off 
size in the continuum description. The effective equa- 
tions of motion are 



dt 
(IX, 
~~dl 
de 
dt 



1 - Pru 



q/ ■ 

— [ei_i(ej_i - o) + ei(e t - a)] - 2/36 

4 in - 6 



Aj Aj + i \ a 

r)i - 7] i+ i - 6 J 2 
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where 8 = 1/Nq. First we describe the gain and loss 
terms in the dynamics of rji . Axons born inside the basin 
of a fascicle contribute to the increase in rji, hence the 
term ftXi (we used Nq = a/ ft to express axon birth rate 
a in terms of ft). A fascicle can lose axons only by in- 
dividual axon deaths, thus the loss term ftrji. Any axon 
which is born in the inter- fascicle empty spaces 6j_i and 
ends up in either of the two neighboring fascicles with 
probability 1/2, hence the gain term (1 /2)/?(ei_i + e.;) 

Next we consider the dynamics of basin size A; . A new 
axon can be born in the gap with a rate a{ti — a) and 
attach to the i-th fascicle with probability 1/2. If it at- 
taches it contributes half the gap size ej/2 towards the 
basin size A;. Hence the gain term (a/A)ei(ei — a). A 
similar contribution to the gain in the basin size comes 
from the other neighboring gap e,_i. The death of a 
boundary axon reduces the basin size by an amount 
SXi/{r]i — 8) (assuming no double occupancy at a lat- 
tice point in the y = level). The contributions of this 
loss coming from two boundaries add up in the total loss 
term 2 x ft8Xi/(r) i — 8). 

Finally, we consider the dynamics of the inter-basin 
gaps £j. The death of boundary axons of neighboring fas- 
cicles i and i+ 1 that border the i-th gap contributes to 
the gain in the gap size. Thus the gain terms ftS\i/(j]i—8) 
and ft5\i+i/(r)i + i — 8). Birth of an axon in the gap re- 
duces the gap size under consideration. The rate of such 
an axon birth is a(ei — a) and on average this event re- 
duces the gap size by an amount ej/2. Thus the loss term 
(a/2)e i (e l - a). 

We use a periodic boundary condition, such that the 
last fascicle is a nearest neighbor of the first fascicle. 
These equations obey the constraint of overall constant 
size • (Xj +e») = 1. It is important to note that the cut- 
off a can be taken to zero meaningfully only after solving 
the differential equations. 

For the simplest non-trivial case involving two fas- 
cicles, the steady state that follows from these equa- 
tions is characterized by r\\ — r\ 2 — 1/2, Ai = A2 = 
1/2 - 5/2 - 26 - dp/16 and e x = e 2 = a/2 + 26 + dp/16. 
We perform a normal mode analysis for small deviations 
from this steady state. The constraint J2i=i 2 (Xi+£i) = 1 
implies that there are only five independent deviations, 
6X1, 6e±, 6e2, 6rji and 8r\ 2 . The linear stability analysis 
about the steady state shows that all the five possible 
modes are stable. Among them, four modes are short- 
lived. For them the deviations decay extremely fast with 
rates ~ ft. However, the fifth mode which in the large- 
size limit can be written as 

(6m = -1, 6D 1 w -1/p, 8E 1 = 0, 8n 2 = 1, 6D 2 w 1/p) 

takes a long time to decay (see Fig. [5]). It involves the 
loss of a boundary axon of one fascicle, which shrinks its 
basin size by 8D\ m —1/p, and simultaneously a gain of a 
boundary axon for the other fascicle, increasing its basin 
size by the equal and opposite amount 6D2 ~ 1/p. This 
operation leaves the inter-basin gap unchanged (5E\ = 0) 



and can be viewed as an exchange of basin space (Fig. [8]). 
The deviations from steady state in this mode decay over 
a very long time scale r ~ n 2 /3/3 where n is the steady 
state value of the fascicle size (n\ — n 2 = n). 

Using the approximate growth of mean fascicle size n ~ 
y b with b = 1/2, we find a power law growth of this time 
scale t ~ y 2b . Notice that the measured time scales r Qp 
and t c obtained from MC simulations show an increase 
with y which approximately obeys the power law y 2b with 
b = 1/2 (see Figf4]). At lower densities, the simulated 
data agrees better with the y 2b power law. In the analytic 
calculation above, we assumed single occupancy of the 
boundary sites of a fascicle basin (removal of a boundary 
axon was assumed to reduce the basin size). At higher 
densities this assumption dose not hold, the boundary of 
a basin does get multiply occupied by axons and thus we 
see a departure from the y 2b power law. 



VI. SINGLE TYPE OF AXONS, WITH 
DETACHMENT 

In the previous section, we concentrated on the "zero- 
temperature" , energy-minimizing dynamics, in which ax- 
ons cannot detach from a fascicle once they become 
part of it. Now we extend our analysis to "finite- 
temperature" Monte-Carlo dynamics, in which the de- 
tachment of growth cones from fascicles become possible. 
It is important to analyze this general case as in the ex- 
perimental studies of fasciculation dynamics [2l|, [2i| [3(| , 
defasciculation events arc clearly observed. An additional 
reason is that the detachment from fascicles is crucial for 
the formation of pure fascicles in a system containing 
multiple axon types - see Section [Villi 

As described in Sec. IIIII in our Monte Carlo sim- 
ulations (with effective temperature set to unity), 
a randomly attempted move to left (right) is ac- 
cepted with probability pl = min[l, exp(— 6E1)] (pn = 
min[l, exp(— 6E r )]) where 6E1 and 8E r are evaluated 
based on the axon occupancy numbers and on the addi- 
tive interaction energy per axon Eh (< 0). The detach- 
ment rate of a growth-cone following a fascicle of size n 
is exp(nEh). Note that for weaker Eh, the detachment 
rate is larger. 



A. Impact of detachment on and r ap 

The approach-to-steady-state time scale r ap decays 
with increasing detachment rates (decreasing inter-axon 
attraction \Eh\). We perform the approach to steady 
state data analysis in the same manner as we did for 
the purely energy-minimizing dynamics discussed in the 
previous section. This analysis gives an estimate of the 
time-asymptotic mean fascicle size rioa as well as the 
approach-to-steady-state time scale T ap as a function of 
y. We perform this analysis for a system of Nq = 100 and 
L = 200 at various strengths of inter-axon interaction Eh- 
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FIG. 9: (Color online) The time-asymptotic fascicle size and 
approach-to-steady-state time scale as a function of y for a 
system with density p = 1/2 (7V = 100 and L = 200) at 
different inter-axon attractions Eh- The time series were 
collected over t = 500T and averaged over 10 3 realizations 
(a) The subtracted time-asymptotic fascicle size n 3 — 
(jioo — c) exp(/3y)/2p c ff follows a power law y 1 ^ 2 . The effective 
density p c g and offset c are treated as fitting parameters with 
c = 5.44±0.62, p cB = 0.52±0.02 for E h =-4,c= 3.44±0.26, 
Pea = 0.53 ± 0.01 for E h = -2, and c = 2.26 ± 0.04, 
pea = 0.41 ±0.002 for E h = -1. (b) The approach-to-steady- 
state time scale r ap gets smaller for weaker attractions Eh, 
however, shows the initial power law growth unless Eh > — 1. 
The dotted line shows a power law y 2b with b — 1/2. Fitting 
errors in r a „ are within 5%. 



Strictly speaking, the presence of detachment invalidates 
our earlier mean field argument for the growth of fascicle 
size with y. While the concept of the basin of a fascicle is 
still meaningful, the basin can "leak", i.e., newly grow- 
ing axons can escape from it through the detachment 
process. Despite this, for most interaction strengths our 
numerical results obey = c+2p c gy 1 / 2 exp(—/3y) with 
the fitting parameter p c g taking the place of p [Fig.[5fa)]. 
The value of p e g is the smallest for the weakest attrac- 
tive interaction plotted in Fig. Uta). For extremely weak 
attraction (e.g., Eh = —0.1) almost no fasciculation can 
occur (p e ff ~ 0) and therefore remains independent 
of y (data not shown) . 

As shown in Fig. \§lb), the time scale T ap for the ap- 



proach to steady state decays with reduced inter-axon 
attraction (i.e., with increased detachment rate). When 
Eh < —2, the growth of r ap with y is similar to the one 
we observed for the "zero-temperature" dynamics. At 
Eh = — 1, however, we find that r ap becomes comparable 
to the mean axonal lifetime T, even though at this value 
of Eh the power-law growth of the mean fascicle size with 
y is still maintained. 

The reduction of r ap with increased detachment rate 
may be understood as a result of enhanced interac- 
tion between neighboring fascicles. As we noted be- 
fore, the slowest time scale for the case of purely energy- 
minimizing dynamics is due to the very slow process of 
exchange of basin size between neighboring fascicles. In 
presence of finite detachment rates, axons from one fas- 
cicle can detach and connect to a neighboring fascicle, 
thereby opening up a new and faster mode of interaction 
between neighboring fascicles. 

B. Impact of detachment on the fascicle size 
distribution 

In this subsection we show that even in presence of 
an appreciable detachment rate (Eh < —1), the scaling 
of fascicle size distribution persists and the scaling func- 
tion retains its overall functional form. However, with 
increasing detachment rates (decreasing l-E^Q, the small 
n portion of the fascicle size distribution gains at the cost 
of bigger fascicles and deviates from the scaling form. 

Fig. QUI shows the steady-state fascicle size distribu- 
tions for a system of L = 800 and Nq = 100. As the 
attraction is decreased from Eh = —4 to Eh = — 1, we 
observe an overall increase in the small n portion of the 
distribution. Except for the lowest values of n, we obtain 
data collapse implying the scaling form 

P s (n,y) = (n)- r cf>(n/(n)) 

where r = 2.1 for E h = -4, -2, and r = 2 for E h = -1. 
Similarly to the case of strictly energy-minimizing dy- 
namics, this scaling is observed only in the intermediate 
range of y values, 10 3 < y < 10 4 , and the scaling function 
is of the form 4>(u) — J\fu exp(— vu~ Xu 2 ) with u = n/ (n). 

The small-n part of the distribution does not scale. 
P s (n,y) is large at n = 1, and drops to lower values 
with increasing n, before it increases again to follow the 
scaling function. The functional form of the initial decay 
of P s (n, y) with n depends on y and also Eh, as we show 
in detail in Fig. HPf rf). 

The discussion in this section shows that there are pa- 
rameter regimes, e.g., the Eh = ~ 1 case discussed above, 
where the emergent time scales are comparable to the 
mean lifetime T of single axons, and at the same time 
the steady state statistics (the mean fascicle size and 
the fascicle size distribution) obey the overall features 
demonstrated by the energy-minimizing dynamics. This 
parameter regime might be utilized to attain a fascicula- 
tion pattern of this type in a relatively short time. 
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FIG. 10: (Color online) Steady state distribution of fascicle sizes P s (n,y) (averaged over 10 4 realizations and the time interval 
10T < t < 25T) for the TVo = 100, L = 800 system at different strengths of inter-axon interactions (a) Eh = —4, (b) Eh = —2, 
and (c) Eh — — 1. The insets show the data collapse. A scaling with B = l/(n) and A = (n) r collapses data obtained for 
y = 1585, 1995, 3162, 5012, 6310, 7943, 10 4 onto a single curve 4>{u) = Nuexp{-vu - Ait 2 ) with u = n/(n) and (a) Af = 262.6, 
v = 0.37, A = 0.74, r = 2.1 (6) M = 318, v = 0.67, A = 0.64, r = 2.1 and (c) JV = 244.5, v = 0.18, A = 1.05, r = 2. The fitting 
to obtain the scaling functions is done above u = 0.4. (d) The same distributions as in (c) presented in a log-log plot. P 3 (n, y) 
decreases with n in the small n regime. The thick line (red) denotes a power law n _5/2 . The arrow denotes the direction 
of increasing y (y = 1, . . . , 10 4 ). Note that the curve obtained at y — 10 4 is the closest to the n _5//2 line at small n. For 
y > 500, the initial decrease of P s (n, y) with n is followed by a subsequent increase which merges with the scaling function. At 
y < 20, the tail of the distribution is exponential (non-interacting limit). Inset of (d): The solid (red) lines are the fascicle size 
distributions for Eh = —0.1 at all j/-levels (10 < y < 10 4 ). It shows a clear single exponential decay characterizing the (almost) 
non-interacting axons. This should be compared with the dotted (blue) lines that show approximate exponential decay of the 
distributions in the range y < 20 for Eh = — 1. 



VII. RELATION TO PARTICLE 
AGGREGATION AND COALESCENCE IN ONE 
DIMENSION 

As we pointed out in Sec. lIVi in a steady-state configu- 
ration at fixed time t, the axon fasciculation with increas- 
ing y may be formally viewed as the evolution of a one- 
dimensional (Id) reaction-diffusion process, where the y 
coordinate takes the meaning of time. As we show in 
this section, this limited analogy can be used to approxi- 
mately understand some steady state properties, e.g., the 
distributions of fascicle sizes and of inter-fascicle spatial 
separations. We stress, however, that the full dynamics 
of our system can not be mapped on to a Id reaction- 
diffusion system. The process of axon turnover, which is 



crucial for the dynamical properties of our system, does 
not have any analog in the Id models we discuss in this 
section. 

A. In absence of detachment 

In this subsection we discuss the relation of our basic 
model, in which axons cannot detach from fascicles, to 
irreversible aggregation and coalescence processes. 

Interpreting the fasciculation of axons with increasing 
y (in a steady state configuration) as Id irreversible ag- 
gregation of particles [mA + nA — > (m + n)A) [l7l . |46I |. 
we find a prediction uexp(— Ait 2 ) (see equation 8.4.24 in 
[17|) for the fascicle-size distribution, which is similar to 
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FIG. 11: (Color online) Histogram of inter-fascicle separation Ax. (a) Rescaling of the histogram H(Ax), obtained at various 
y-levels indicated in the legend, leads to data collapse. The line through the collapsed data is a function VAxexp(— QAx 2 ) with 
V = 0.0085 and Q = 0.0002. (6) The scale factors A, B obtained at different y-levels show power-law dependence A ~ y ' 96 
and B ~ y~ 0A7 indicated by the lines through the data points. 



the true distribution uexp(— vu — Am 2 ) (Sec. IV Bp [TTj j . 
Note that the distribution obtained from the mapping 
to irreversible aggregation lacks the exponential part 
exp(— vu). As we explained in Ref. having v = 
in our model would require r ap — oo. The absence of the 
exponential part is therefore consistent with the absence 
of turnover-based dynamics in the Id analogy. 



The merging of fascicles with growing y may alterna- 
tively be interpreted as an irreversible coalescence process 
A + A — > A [ill, viewing each fascicle as a particle A. 
The pattern formation in Id irreversible coalescence had 
been quantified by the inter-particle distribution func- 
tion (IPDF) [III, the distribution of distance between 
neighboring particles. The steady state of this process 
is trivial, a completely empty space. The IPDF is ob- 
tained at finite time t before this steady state arrives, 
and has the form (x/AT>t) exp(— x 2 /8T)t) where D de- 
notes the particle diffusion constant ■ The change in 
the inter-fascicle separation distribution with increasing 
y in our model may be viewed as equivalent to the time 
evolution of IPDF in irreversible coalescence. This leads 
us to a prediction of the distribution of inter-fascicle sep- 
aration H(Ax,y) s» (Ax/y) exp(— p Ax 2 /y). Using the 
same stochastic simulation that we used to obtain Fig. [51 
we calculated the histogram of spatial separations Ax 
between fascicles identified at various w-levels. Note that 
the separation between two neighboring fascicles Ax is 
measured at the same y-level at which the fascicles are 
identified, and is different from the gap between fasci- 
cle basins E (shown in Fig. [2]). A rescaling of H(Ax) 
by a factor A(y) and Aa; by B(y) leads to a data col- 
lapse (Fig. Illf a)), and approximate power law depen- 
dences A ~ y p and B ~ y~ q with p — 0.96 and q = 0.47 
(Fig. fTTT b)). This result is in reasonable agreement with 
the above-mentioned form of H(Ax,y) that predicts a 
scaling function? 3 Ax exp(—QAx 2 ) (Fig.fTTTa)). and scal- 
ing exponents p — 1 and q — 1/2. 



B. In presence of detachment 

In this subsection we discuss analogies of our model of 
axon fasciculation in presence of detachment to reversible 
aggregation and coalescence processes in Id. Due to spe- 
cific features of our model, only qualitative analogies to 
models from the Id literature can be made. 

Detachment events are partially captured when the 
axon fasciculation with increasing y in a fixed-time con- 
figuration is formally viewed as a Id diffusion with re- 
versible aggregation, the chipping model: mA + nA — > 
(to + n)A and mA — > (m — 1)A + A [la| . denoting each 
axon by a particle A. (Note, however, that steady state 
configurations of our model show splitting of fascicles 
with increasing y (e.g., see Fig. [DJb)) into two fasci- 
cles containing multiple axons. The chipping model does 
not include the analog of such a process). The chipping 
model posseses a non-trivial steady state (t — > oo, cor- 
responding to y — > oo within our model). It shows a 
dynamic phase transition associated with particle den- 
sity p O [l!|. The steady state distribution of clus- 
ters of size n is predicted to be P s {n) ~ exp(— n/n*) at 
p < p c — VI + w — 1, where w denotes a constant single 
particle chipping rate. At the critical density p = p c , 
the distribution changes its shape to P s (n) ~ n™ 5 / 2 . At 
density above p c this power-law distribution remains un- 
altered, and in addition to the power law distributed 
clusters one gets a single cluster of diverging size [l8| . 
Fig. [10] shows fascicle size distributions obtained from our 
model at various interaction strengths Eh- Note that, in 
contrast to the chipping model, in our model the rate 
w = 1/[1 + cxp(— nEh)\, with which a single axon de- 
taches from a fascicle, depends on the fascicle size n. 
Only for the weakest interaction Eh = —0.1 (inset of 
Fig. [TOT d)). the detachment rate from a two-axon fasci- 
cle w = 0.45 corresponds to a critical density p c = 0.2 
which is greater than the axon density p = 1/8. The cor- 
responding fascicle size distribution shows a form con- 



13 



sistent with ~ exp(— n/n*) (inset of Fig. [TO^d)). For 
Eh = — 1 the detachment rate from a two-axon fascicle 
w = 0.12 corresponds to a critical density p c = 0.06 < p 
(= 1/8). The fascicle size distribution in the region of 
small n, for Eh = — 1, shows rough agreement with the 
power law ~ n~ 5 / 2 (Fig. noT d)). This change in shape of 
the fascicle size distribution from an exponential decay 
to a power-law decay at small n, thus, is consistent with 
the dynamical phase transition predicted by the chipping 
model. Since in our model the detachment rate of axons 
exponentially decays with fascicle size n, for larger fasci- 
cles the detachment rate w gets so small that the analogy 
with the reversible chipping model breaks down, and the 
behavior of the system becomes analogous to irreversible 
aggregation. The fascicle size distribution in the region 
of larger n fFig.UOp becomes indistinguishable from axon 
fasciculation in absence of detachment. 

We note that, for our model in presence of detachment, 
the change in the distribution of inter-fascicle spatial sep- 
aration with increasing y can not be easily understood 
in terms of reversible coalescence A + A =; A [l6l . l48j 
(denoting each fascicle as a particle A). The main rea- 
sons are: (i) the reverse reaction A — > A + A allows 
for splitting of a single axon into two, a mechanism not 
allowed in our model; (ii) in contrast to reversible coa- 
lescence, the probability of splitting of a fascicle in our 
model decays rapidly with increasing fascicle size. Thus 
the y-independent distribution of inter-fascicle separation 
c s exp(— c s Aa;), expected from IPDF of reversible coales- 
cence [Ty] , is never reached. At large y, we find a dis- 
tribution of inter-fascicle separation that conforms more 
to VAx exp(— Q Ax 2 ) (data not shown), consistent with 
irreversible coalescence (see previous subsection). 

We note again that time-dependent quantities in our 
model have no analog in the mapping to the Id models 
we discussed above. The emergence of density-de pen dent 
long time scales in reversible coalescence [HI, |48| . l49j 
therefore has no relation to the long time scales in our 
model, which are due to a very slow reorganization of fas- 
cicle basins (Sec. lVD|) — a consequence of axon turnover. 



VIII. MIXED POPULATION OF MULTIPLE 
AXON TYPES 

The axons of olfactory sensory neurons expressing dis- 
tinct odorant receptors [36f are believed to have short- 
range interactions with interaction strengths that are 
correlated with the type of receptors the neurons ex- 
press [4l|, |42|, [5^, HH- I n the framework of our model, 
this corresponds to the introduction of multiple types of 
random walkers, with type-dependent probabilities for 
attachment to / detachment from fascicles. Neuronal 
systems containing multiple types of axons are known to 
achieve pure and stable connections. In the olfactory sys- 
tem, it may be expected [HJ] that when the growth cone 
is located in a relatively pure environment (i.e., when it is 
in contact with axons mainly of its own type) , it leads to 



reduced axonal turnover. While we do not include such 
effects in our model, in this section we evaluate the mean 
purity of axon environment S to characterize the sorting 
dynamics. 



A. Mean purity 

Let us consider a system containing two types of axons 
named r and b. Assume that a fascicle contains n r of r- 
axons and n b of b- axons. If the i-th axon in the fascicle 
is of type r, the purity of environment that this axon 
encounters within the fascicle is Si = (n r — n b )/ (n r + n b ), 
while if the axon is of type 6, the purity of environment 
is Si = (n b — n r )/(n r + n b ). Then the mean purity of 
environment obtained by averaging over all N(y) axons is 

S = (X}i=i S i)/N(y)- The partial sum within a fascicle 
gives n r (n r - n b )/(n r + n b ) + n b {n b - n r )/(n r + n b ) = 
(n r — n b ) 2 /(n r + n b ). Thus, the degree of sorting can be 
quantified as the mean purity of axon environment 
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Notice that < S < 1; S = 1 corresponds to completely 
pure fascicles containing only one type of axons, whereas 
5 = describes fascicles containing an equal mixture of 
the two axon types. 



B. Approach to steady state 

We consider a system of size L — 800 having TVq = 50 
axons of type r and Nq = 50 axons of type b at the y = 
level (p = (N£ + N$)/L = 1/8). We use the homotypic 
interaction energy (interaction between r-r or b-b) Eh 

- and the heterotypic interaction energy (between the 
two different types r-b) E a = —0.1. We monitor the time 
evolution of the mean purity of environment S and the 
mean fascicle size n. A typical configuration is shown in 
Fig. W[b). The mean fascicle size reaches the steady state 
value at each y-level with two characteristic time scales, 
similarly to the case of a system containing only one type 
of axons (Fig. [3]). The shorter time scale is intrinsic, 
the mean lifetime of a single axon T, and the other one 
is the emergent approach-to-steady-state time scale r n 
(equivalent to r ap defined in the caption of Fig. [3]). As in 
the case of a system containing only one type of axons, 
this time scale r„ can be orders of magnitude larger than 
7 (Kig.[I3-/::. 

The intrinsic time scale T does not appear in the dy- 
namics of S, however. S approaches its steady state 
value 6*00 with a single time scale t s . The measured 
mean purity S fits to the form S — Soo + rexp(— t/r s ) 
(Fig. HHa)). The steady-state mean purity initially 
grows with y, reaching a maximum beyond which 
decreases (Fig. [T27 M). This non-monotonic behavior is 
seen both for low (p = 1/8) and high (p — 1/2) density 
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FIG. 12: (Color online) This plot shows the approach to steady state for mean fascicle size and mean purity of fascicles in 
a system containing iVJ = 50 r-axons and N$ = 50 6-axons that interact through homotypic interaction strength Eh = —4 
and heterotypic interaction strength E — —0.1. Two mean axon densities p = 1/8 (system size L — 800) and p = 1/2 
(L = 200) are used, (a) Approach-to-steady-state data for mean purity 5 as a function of time t/T in a system with density 
p — 1/8. All the data were collected over 500T and averaged over 10 3 realizations . The data at each y-level fits to the form 
S(t) — Soo + r exp(— t/r s ) where 5*00 is the asymptotic mean purity and t s is the time scale of approach to steady state. The 
fitting is shown for the data set at y = 10 4 , where Soo = 0.73, r = 0.19 and t s = 63.64 with all the fitting errors being less 
than 2%. (b) Asymptotic mean purity Soo as a function of y at p — 1/8 and p = 1/2. Fitting errors in Soo are within 3%. 
(c) Asymptotic subtracted mean fascicle size n s = (rioo — c) exp(/3y) /2p c g as a function of y. This follows y b with b = 1/2. 
p eS is treated as a fitting parameter. For p = 1/2, c = 1.88 ± 0.2 and p cB = 0.392 ± 0.005. For p = 1/8, c = 1.13 ± 0.05 
and p c ff = 0.059 ± 0.001. (d) Approach-to-steady-state time scales, r s for mean purity S and r n for mean fascicle size n, as 
a function of y for systems with p — 1/8 and p = 1/2. Fitting errors in r s and t„ are within 3%. The time scales show an 
approximate power law growth y 2b with 6=1/2 denoted by the solid black line. Data are shown up to a t/-level where the time 
scales extracted from the fitting procedures remain less than half the total run time t — 500T. 



systems (Fig. [T2T &)). Thus for the combination of strong 
homotypic and weak heterotypic attractive interaction, 
the fascicles achieve the highest purity at a particular dis- 
tance from their starting point, beyond which the typical 
fascicle keeps on losing purity. From energetic considera- 
tions, one may understand this behavior as follows. The 
high strength of the homotypic attraction compared to 
the heterotypic one leads to sorting and thus the ini- 
tial growth in purity at lower y-levels. Once one obtains 
highly pure and sufficiently large fascicles, however, the 
heterotypic interaction will merge the r-dominated and 
fo-dominated fascicles to form larger mixed fascicles, and 
thereby lead to a loss in mean purity. The effect of the 
weak heterotypic interaction becomes significant only at 



higher y-levels, where large fascicles are formed (recall 
that in our model, the interactions are additive). 

Similar to the case of a system containing only a sin- 
gle axon type in presence of detachment, we find that 
the steady state mean fascicle size grows as rioo = c + 
2p e ff2/ 1 / 2 exp(— j3y). Fig. [T2f c) shows that the subtracted 
mean fascicle size n s = (rioo — c) exp(/3y)/2,o c ff follows the 
power law y b with b = 1/2. The effective density p e s and 
offset fascicle size c are treated as fitting parameters. 

The time scales for approach to steady state r„ and 
t s grow with y following an approximate power law y 2b 
with 6 = 1/2 [Fig. IT27 <i)]. This behavior is seen to be 
independent of density p, in contrast to the single-type 
case, where we found reliable y 2b growth of the emergent 
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FIG. 13: (Color online) The distributions of fascicle composition P(n r ,n b ) calculated by averaging over 10 4 realizations and 
time interval 10T < t < 15T in a system with homotypic interaction Eh = —4, heterotypic interaction E = —0.1, number of 
axons Nq = N$ = 50 and system size L = 200. The plots show the distributions calculated at y = 10, 10 2 , 10 3 , 10 4 . 



time-scales only at lower densities. Recall that, reduced 
inter-axon interaction strengths lead to lower effective 
densities p e g (Fig. Oa)). In the present case of mixed 
population of axons of two types, the heterotypic inter- 
action is very weak and may have lead to the effectively 
low-density (y 2b ) power-law growth. We note that r n and 
t s turn out to be approximately equal to each other at 
all y-levcls. 

C. Distribution of fascicle composition 

In this section we briefly discuss the distribution of 
fascicle composition P(n r ,n b ), measured as the number 
of fascicles with n r r-axons and n b 6-axons at a speci- 
fied y-level and time t. In Fig. [T3] we show this distri- 
bution obtained by collecting data within the time in- 
terval 10T < t < 15T at various y-levels for a system 
of N£ = N b = 50 and L = 200 (p = 1/2). From the 
approach-to-steady-state data (Fig. Q2]) it is clear that 
at t = 10T, steady state is reached only up to y ss 100. 
Thus the distributions obtained at y = 10 3 , 10 4 in Fig. [13] 
are far from steady state [55| . All the plots in Fig. [T51 
show a pronounced maximum for evenly mixed fascicles, 
meaning that most of the fascicles we obtain are mixed by 



type. However, a careful look at the plots reveals off-peak 
features of the distribution with reasonable weight that 
reflect the presence of fascicles with highly asymmetrical 
composition (e.g., along the n b = 1 line for the plot at 
y = 100 in Fig. 1X31 the maximum of the distribution is 
at n r — 5). 

Finally, we comment on repulsive vs. attractive het- 
erotypic interactions. The attractive heterotypic interac- 
tion generates larger fascicles asymptotically, but induces 
a reduction of purity at large y. A repulsive heterotypic 
interaction (combined with attractive homotypic interac- 
tion) would generate enhanced sorting, however it would 
be at the cost of a decreased mean fascicle size at all y 
levels. 



IX. SUMMARY 

In this paper, we provided a simple model of the dy- 
namics of axon fasciculation and sorting. To allow us 
to concentrate on the collective effects that arise from 
axon-axon interactions in a large population of axons, we 
chose a particularly simple implementation of single axon 
growth. In our model, each growing axon is represented 
as a directed random walk (Sec. IIII A|) . The common 
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preferred growth direction may arise, e.g., from the in- 
fluence of a spatially distributed guidance cue emitted by 
a distant target. Other than this common influence, we 
do not include in our model the guidance of axon growth 
cones by graded guidance cues, and restrict our attention 
to axon-axon interactions. The interaction of a growth 
cone with other axons is modeled as a short-range attrac- 
tive interaction between a random walker and the trails 
of other random walkers (Sec. IIII A"l) . In addition, we in- 
corporated neuronal turnover (characteristic of, e.g., the 
mammalian olfactory system) by assigning a finite life- 
time to each growing axon (Sec. IIII Bj) . 

The strength of the axon-axon interaction was 
parametrized as an effective energy in a Monte Carlo up- 
date. In the energy-minimizing dynamics (corresponding 
to zero effective temperature in the Monte Carlo scheme), 
once a growth cone attaches to an axon fascicle, it will 
never detach (become a free random walker) again. For 
such dynamics, we extended our previous numerical and 
analytical results of Ref. [ll[ (Sec. [Vj . In the general dy- 
namics (corresponding to unit effective temperature), ax- 
ons may detach from fascicles, with a rate that increases 
with decreasing axon-axon interaction strength. We sys- 
tematically studied how such detachment events modify 
the basic dynamics with no detachments (Sec. I VI I) . Fi- 
nally, we investigated a system with two types of axons, 
and analyzed the sorting dynamics arising from a strong 
homotypic interaction (between axons of same type) com- 
bined with a weak heterotypic interaction (between axons 
of different types) (Sec. IVIII[) . Our main findings were 
as follows. 

The tendency to fasciculate is reflected in the growth 
of the mean fascicle size n with the distance y in the 
preferred growth direction of the axons. Using a mean 
field argument we showed that for the energy-minimizing 
dynamics the mean fascicle size n should grow as n ~ 
2y 1 / 2 pexp(-/3y). This agrees with the numerical results 
of Sec. El In Sec. Ell and IVTiTl we showed that this 
growth law persists even in presence of detachment, with 
the average axon density p replaced by a fitting parame- 
ter p cS . 

A more detailed characterization of the steady state 
is a position-dependent fascicle size distribution. Within 
the scaling regime L <C y -C (L/2) 2 , this distribution 
obeys a scaling law P s (n,y) = (n(y))~ r (f>(n/ (n(y))) with 
r = 2.1 and the scaling function 4>{u) = ACwexp(— vu — 
Xu 2 ) (SecEJ- At higher y-levels the distribution becomes 
bimodal, a new maximum arises which is characteristic 
of the complete fasciculation. Even in the presence of de- 
tachment, the scaling behavior of fascicle size distribution 
remains valid over a wide range of interaction strengths 

(Sec.EUandEEEIl)- 

The dynamics of reorganization of fascicles at high y- 
levels was found to be extremely slow. The emergent time 
scales, e.g., the approach-to-steady-state time r ap or the 
auto-correlation time at steady state r c can be orders 
of magnitude larger than the mean lifetime of an axon 
T. In Sec. IV Dl using an analytical model of effective 



dynamics involving two neighboring fascicles, we showed 
that the slowest mode of this dynamics corresponds to 
the exchange of basin space between the two fascicles, 
and grows with distance as t ~ y. This behavior of time 
scales survives even in the presence of detachment (shown 
in Sec. lVTAl . 

While our model is two-dimensional, some limited 
analogies can be made to one dimensional (Id) models 
of aggregation, coalescence, and chipping (Sec. VII). We 
introduced a mapping in which the progressive fascicu- 
lation with increasing y (at a fixed time) in our model 
is mapped on to the time evolution within a Id sys- 
tem of interacting particles. Using results from the lit- 
erature on Id models, we then obtained predictions for 
stationary quantities in our model. Thus interpreting 
each axon as a particle A in the irreversible aggrega- 
tion model mA + nA —> (m + n)A [13, EH, we obtained 
the prediction uexp(— Xu 2 ) for the distribution of fasci- 
cle sizes, which is similar to the true steady-state dis- 
tribution wexp(— vu — Xu 2 ) in our model with energy- 
minimizing dynamics. Likewise, interpreting each fas- 
cicle of axons as a particle A in the irreversible coales- 
cence model A + A — > A [l6| , we obtained the predic- 
tion (Ax/y) exp(— uAx 2 /y) for the distribution of sep- 
arations between fascicles, which agrees approximately 
with numerical results from our model (Sec. IVIlT) . A lim- 
ited analo gy can also be made between the Id chipp ing 
(Ref. [H, H3) or reversible coalescence (Ref. [l6l l48j) 
models and our model in the presence of detachment. 
Since in our model, the rate of detachment decreases with 
the fascicle size, the reversible interaction models are rel- 
evant only at low axon-axon interaction strengths. In this 
range of parameters, we were able to relate the observed 
changes of distribution of fascicle sizes in our model to 
the phase transition that occurs in the chipping model of 
Ref. [IH (see Sec. IVIII) . We stress again, however, that 
the mapping to these Id models can say nothing about 
the time-dependent quantities in our model (time in our 
model has no analog in the Id models). Therefore, for 
example, the slow time scales discussed in the reversible 
coalescence model of Refs. [ljl Eg], Ei| are unrelated to 
the slow time scales present in our model. 

In Sec. lVIlIl we analyzed a system with two types of ax- 
ons and type-specific interactions. In this system, axons 
sort into fascicles according to axon type. We quanti- 
fied the degree of sorting by introducing the mean pu- 
rity S of axon environment within fascicles. For the case 
of strongly attractive homotypic interaction and weakly 
attractive heterotypic interaction, we showed that the 
degree of sorting S varies with distance y in a non- 
monotonic manner and has a single maximum. 



X. OUTLOOK 

In this paper we analyzed a model aiming to describe 
the formation of axon fascicles and the sorting of fasci- 
cles by neuronal types in the mammalian olfactory sys- 
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tern. Our goal was to systematically investigate the gen- 
eral non-equilibrium statistical mechanics aspects of the 
model, leaving the task of building quantitative connec- 
tions with physiology for the future. 

To conclude, we discuss possible generalizations of the 
basic model defined in this article, and the applicabil- 
ity to biological data on axon fasciculation. First, we 
note that in our discussion of the dynamical properties 
of the system, it was essential that the random walkers 
moved in two (rather than three) spatial dimensions, and 
therefore cannot cross each other without interacting. In 
contrast, in a three dimensional system, the concept of 
fascicle basins would lose its validity. The resulting fas- 
cicle dynamics in three dimensions is expected to be sig- 
nificantly different from the two-dimensional dynamics, 
which we showed to be governed by the competition of 
fascicles for basin space. It will be necessary to examine 
to what extent the assumption of two dimensionality is 
satisfied in the olfactory system. However, this assump- 
tion is effectively satisfied in studies of growth in neu- 
ronal cell culture [U |43|, [52j , in which the axons move 
on a plane surface, and interact when crossing each other. 
In Ref . [I3J , a fluorescence- based method is proposed for 
extracting the distribution of fascicle sizes; such experi- 
ments would permit a direct test of our model. Note that 
in some of these studies, our assumption of very strong 
adhesion of axons to the substrate is not satisfied. Events 
not included in our model, such as the gradual straight- 
ening of axon shafts or the local zippering / unzipper- 
ing of fasciculated axons (as observed in Ref. [52[) may 
therefore occur. Note also that natural boundary con- 
ditions in such cell cultures are either free boundary [43[ 
or confining channels [U [52j , in contrast to the periodic 
boundary condition used in our simulations. 

Recall that the mouse olfactory system contains about 
1000 axon types. In this paper, we considered only up to 
2 types of axons. It would be interesting to examine if 
qualitatively new features emerge in systems with many 
axon types and a range of heterotypic interactions. This 
would require, however, significantly longer simulations. 



Recently, Ref. |13[ examined the role of axon-axon 
interactions in the fasciculation and sorting of axons 
belonging to mouse olfactory sensory neurons. This 
showed that inter-axon repulsive interactions arising from 
Neuropilin-Scmaphorin signaling play an important role 
in the axon sorting. In the pre-target region (before 
reaching the olfactory bulb) , the amount of sorting grows 
with distance from olfactory epithelium [131 ]. The pre- 
target axon sorting is shown to affect the topographic 
map formation by the neurons in the olfactory bulb [13l — 
UU. Using a mutant mouse, Ref. [HI further showed 
that heterotypic axons sort even in absence of the ol- 
factory bulb, i.e., in complete absence of axon-target in- 
teractions. This kind of experiments forms a suitable 
ground for the application of our model. Notice that in 
our model, we have shown that a weaker heterotypic at- 
traction combined with a stronger homotypic attraction 
already leads to sorting. An effective repulsion between 
the two types would enhance the amount of sorting, how- 
ever, it would be at the expense of the size of the fascicles 
formed. 

Our immediate future goal is to extract the parameter 
values of our model from controlled in vitro experiments 
and then use the model (and its possible extensions) to 
analyze in vivo data on olfactory pattern formation in 
mice. 
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